This is the pipeline used to analyze the irCLIP-RNP TMT datasets for 13 RBPs from three different gel sections ranging from 60-120kDa, 120-225kDa, and 225-350kDa. The experiment was performed in HEK293T and HepG2 cells.

1. Prepare the dataset

#Needed libraries
library(DEP2)
library(tidyverse)
library(ggplot2)
library(data.table)
library(pheatmap)
library(RColorBrewer)
library(gplots)
library(hrbrthemes)
library(pacman)
library(textshape)
library(ggExtra)
library(viridis)
library(purrr)
library(hexbin)
library(DESeq2)
library(ggpubr)
library(UpSetR)
library(dplyr)
library(Clipper)
library(factoextra)
library(paletteer)
library(corrplot)
library(psych)
library(ggpmisc)
library(gprofiler2)
library(viridis)
library(GGally)
library(igraph)
library(rstatix)
library(limma)
# Open TMT data that were searched with MaxQuant and processed with Perseus
ABCF1_data <- read.delim("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/ABCF1_BZ59.txt")
DDX5_data <- read.delim("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/DDX5_BZ6.txt")
FUS_data <- read.delim("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/FUS_BZ3.txt")
hnA2B1_data <- read.delim("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/hnA2B1_BZ55.txt")
hnC_data <- read.delim("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/hnC_BZ56.txt")
hnM_data <- read.delim("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/hnM_BZ57.txt")
hnU_data <- read.delim("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/hnU_BZ58.txt")
ILF2_data <- read.delim("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/ILF2_BZ1.txt")
ILF3_data <- read.delim("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/ILF3_BZ2.txt")
NAT10_data <- read.delim("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/NAT10_BZ60.txt")
NONO_data <- read.delim("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/NONO_BZ4.txt")
RBFOX2_data <- read.delim("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/RBFOX2_BZ86.txt")
SFPQ_data <- read.delim("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/SFPQ_BZ5.txt")

#get the unique gene names and protein IDs
colnames <- c("name", "ID", "H293T.noUV.R1.L",  "H293T.noUV.R1.HM", "H293T.UVC.R1.L",   "H293T.UVC.R1.M",   "H293T.UVC.R1.H",   "H293T.UVC.R2.L",   "H293T.UVC.R2.M",
              "H293T.UVC.R2.H", "HepG2.noUV.R1.L",  "HepG2.noUV.R1.HM", "HepG2.UVC.R1.L",   "HepG2.UVC.R1.M",   "HepG2.UVC.R1.H",   "HepG2.UVC.R2.L",   
              "HepG2.UVC.R2.M", "HepG2.UVC.R2.H")
data_prep <- function(data) {
  data$Gene.names <- str_match_all(data$Fasta.headers, "GN=(.*?) PE") %>%  lapply(., function (x) str_c(x[,2],collapse='; ')) %>%  unlist()
  data$Prot.IDs <- str_match_all(data$Fasta.headers, "(?<=sp\\|)[[:alnum:]]+") %>%  lapply(., function (x) str_c(x[,1],collapse='; ')) %>%  unlist()
  data$Gene.names %>% duplicated() %>% any() # check for duplicates
  data <- make_unique(data, "Gene.names", "Prot.IDs", delim = ";")
  data$name %>% duplicated() %>% any() # must be false
  
  data_f <- data[,c(tail(grep("name", colnames(data) ), 1),tail(grep("ID", colnames(data) ), 1),grep("Reporter", colnames(data) ))]
  data_f$name[data_f$name == "RBFOX1"] <- "RBFOX2"
  data_f$ID[data_f$ID == "Q9NWB1"] <- "O43251"
  return(data_f)
}

#Open TMT data
ABCF1 <- data_prep(ABCF1_data)
colnames(ABCF1) <- colnames
DDX5 <- data_prep(DDX5_data)
colnames(DDX5) <- colnames
FUS <- data_prep(FUS_data)
colnames(FUS) <- colnames
hnA2B1 <- data_prep(hnA2B1_data)
hnA2B1 <- hnA2B1[,c(1:3,7,5:6,4,8:18)]
colnames(hnA2B1) <- colnames
hnC <- data_prep(hnC_data)
colnames(hnC) <- colnames
hnM <- data_prep(hnM_data)
hnM <- hnM[,c(1:7,10,8,9,11:18)]
colnames(hnM) <- colnames
hnU <- data_prep(hnU_data)
colnames(hnU) <- colnames
ILF2 <- data_prep(ILF2_data)
colnames(ILF2) <- colnames
ILF3 <- data_prep(ILF3_data)
colnames(ILF3) <- colnames
NONO <- data_prep(NONO_data)
colnames(NONO) <- colnames
NAT10 <- data_prep(NAT10_data)
NAT10 <- NAT10[,c(1:10,11,16,13,14,12,18,15,17)]
colnames(NAT10) <- colnames
RBFOX2 <- data_prep(RBFOX2_data)
colnames(RBFOX2) <- colnames
SFPQ <- data_prep(SFPQ_data)
SFPQ <- SFPQ[,c(1:15,17,16,18)]
colnames(SFPQ) <- colnames
SFPQ
#Save intensities
write.table(ABCF1, file ="~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/ABCF1_TMT_intensities.txt", row.names = FALSE, sep = '\t', quote = FALSE)
write.table(DDX5, file ="~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/DDX5_TMT_intensities.txt", row.names = FALSE, sep = '\t', quote = FALSE)
write.table(FUS, file ="~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/FUS_TMT_intensities.txt", row.names = FALSE, sep = '\t', quote = FALSE)
write.table(hnA2B1, file ="~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/hnA2B1_TMT_intensities.txt", row.names = FALSE, sep = '\t', quote = FALSE)
write.table(hnC, file ="~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/hnC_TMT_intensities.txt", row.names = FALSE, sep = '\t', quote = FALSE)
write.table(hnM, file ="~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/hnM_TMT_intensities.txt", row.names = FALSE, sep = '\t', quote = FALSE)
write.table(hnU, file ="~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/hnU_TMT_intensities.txt", row.names = FALSE, sep = '\t', quote = FALSE)
write.table(ILF2, file ="~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/ILF2_TMT_intensities.txt", row.names = FALSE, sep = '\t', quote = FALSE)
write.table(ILF3, file ="~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/ILF3_TMT_intensities.txt", row.names = FALSE, sep = '\t', quote = FALSE)
write.table(NAT10, file ="~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/NAT10_TMT_intensities.txt", row.names = FALSE, sep = '\t', quote = FALSE)
write.table(NONO, file ="~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/NONO_TMT_intensities.txt", row.names = FALSE, sep = '\t', quote = FALSE)
write.table(RBFOX2, file ="~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/RBFOX2_TMT_intensities.txt", row.names = FALSE, sep = '\t', quote = FALSE)
write.table(SFPQ, file ="~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/SFPQ_TMT_intensities.txt", row.names = FALSE, sep = '\t', quote = FALSE)

2. DEP analysis for UVC enriched proteins

We used the following design to create a SummarizedExperiment.

# Load design matrix
design <- read.delim("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/0_design.txt")
design

We determined RDAPs by comparing UVC and no-UV samples. Only proteins with FDR < 0.1 and FC > 1.2 were selected.

#Function to perform DEP analysis
DEP_analysis_UVC <- function (data) {
ecols <- c(grep("H293T", colnames(data)), grep("HepG2", colnames(data)) )
data[,3:18] <- 2^data[,3:18]
se <- make_se(data, columns = ecols, expdesign = design)
se <- normalize_vsn(se)
se_H293T <- se[,se$celltype == "H293T"]
diff_H293T <- test_diff(se_H293T, type = "control", control = "noUV_H293T", fdr.type = "BH")
dep_H293T <- add_rejections(diff_H293T, alpha = 0.1, lfc = log2(1.2))
results_H293T <- get_results(dep_H293T)
results_H293T$HEK293T_sign <- ifelse((results_H293T$high_cut_H293T_vs_noUV_H293T_significant == "TRUE" |
                              results_H293T$medium_cut_H293T_vs_noUV_H293T_significant == "TRUE" |
                              results_H293T$low_cut_H293T_vs_noUV_H293T_significant == "TRUE"), "HEK293T", "none")

se_HepG2 <- se[,se$celltype == "HepG2"]
diff_HepG2 <- test_diff(se_HepG2, type = "control", control = "noUV_HepG2", fdr.type = "BH")
dep_HepG2 <- add_rejections(diff_HepG2, alpha = 0.1, lfc = log2(1.2))
results_HepG2 <- get_results(dep_HepG2)
results_HepG2$HepG2_sign <- ifelse((results_HepG2$high_cut_HepG2_vs_noUV_HepG2_significant == "TRUE" |
                              results_HepG2$medium_cut_HepG2_vs_noUV_HepG2_significant == "TRUE" |
                              results_HepG2$low_cut_HepG2_vs_noUV_HepG2_significant == "TRUE"), "HepG2", "none")
results <- merge(results_H293T, results_HepG2[,-2], by = "name")
results$UVC <- ifelse((results$HEK293T_sign == "HEK293T" | results$HepG2_sign == "HepG2") , "UVC", "none")

return(list(se = se, HEK293T_res = results_H293T, HepG2_res = results_HepG2, all = results))
}

ABCF1_DEP_UVC <- DEP_analysis_UVC(ABCF1)
DDX5_DEP_UVC <- DEP_analysis_UVC(DDX5)
FUS_DEP_UVC <- DEP_analysis_UVC(FUS)
hnA2B1_DEP_UVC <- DEP_analysis_UVC(hnA2B1)
hnC_DEP_UVC <- DEP_analysis_UVC(hnC)
hnM_DEP_UVC <- DEP_analysis_UVC(hnM)
hnU_DEP_UVC <- DEP_analysis_UVC(hnU)
ILF2_DEP_UVC <- DEP_analysis_UVC(ILF2)
ILF3_DEP_UVC <- DEP_analysis_UVC(ILF3)
NAT10_DEP_UVC <- DEP_analysis_UVC(NAT10)
NONO_DEP_UVC <- DEP_analysis_UVC(NONO)
RBFOX2_DEP_UVC <- DEP_analysis_UVC(RBFOX2)
SFPQ_DEP_UVC <- DEP_analysis_UVC(SFPQ)
write.table(ABCF1_DEP_UVC$all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/2_UVC_enriched/ABCF1_TMT_DEP_UVC_results.txt", row.names = FALSE, sep = "\t")
write.table(DDX5_DEP_UVC$all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/2_UVC_enriched/DDX5_TMT_DEP_UVC_results.txt", row.names = FALSE, sep = "\t")
write.table(FUS_DEP_UVC$all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/2_UVC_enriched/FUS_TMT_DEP_UVC_results.txt", row.names = FALSE, sep = "\t")
write.table(hnA2B1_DEP_UVC$all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/2_UVC_enriched/hnA2B1_TMT_DEP_UVC_results.txt", row.names = FALSE, sep = "\t")
write.table(hnC_DEP_UVC$all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/2_UVC_enriched/hnC_TMT_DEP_UVC_results.txt", row.names = FALSE, sep = "\t")
write.table(hnM_DEP_UVC$all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/2_UVC_enriched/hnM_TMT_DEP_UVC_results.txt", row.names = FALSE, sep = "\t")
write.table(hnU_DEP_UVC$all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/2_UVC_enriched/hnU_TMT_DEP_UVC_results.txt", row.names = FALSE, sep = "\t")
write.table(ILF2_DEP_UVC$all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/2_UVC_enriched/ILF2_TMT_DEP_UVC_results.txt", row.names = FALSE, sep = "\t")
write.table(ILF3_DEP_UVC$all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/2_UVC_enriched/ILF3_TMT_DEP_UVC_results.txt", row.names = FALSE, sep = "\t")
write.table(NAT10_DEP_UVC$all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/2_UVC_enriched/NAT10_TMT_DEP_UVC_results.txt", row.names = FALSE, sep = "\t")
write.table(NONO_DEP_UVC$all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/2_UVC_enriched/NONO_TMT_DEP_UVC_results.txt", row.names = FALSE, sep = "\t")
write.table(RBFOX2_DEP_UVC$all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/2_UVC_enriched/RBFOX2_TMT_DEP_UVC_results.txt", row.names = FALSE, sep = "\t")
write.table(SFPQ_DEP_UVC$all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/2_UVC_enriched/SFPQ_TMT_DEP_UVC_results.txt", row.names = FALSE, sep = "\t")

#3. Two-factor analysis between cell types and gel sections To analyze the cell-type differences, we used the limma package with a two-factor design (cell type and section).

#Perform DEP analysis using a two-factor model
extract_results_table <- function(fit) {
  tests <- colnames(fit$coefficients) %>%
    set_names(., .)
  map(tests, \(x) {
    tibble(gene = rownames(fit$coefficients),
           logFC = fit$coefficients[,x],
           s2.post = fit$s2.post,
           p.value = fit$p.value[,x]) %>%
      mutate(fdr = p.adjust(p.value, method = 'BH'))
  })
}

DEP_analysis <- function (se, rbp, rbp2) {
se_UVC <- se[,se$crosslinking == "UVC"]

#Get contrasts for each cell line and each normalization
model_vsn <- model.matrix(~ celltype * section, colData(se_UVC))
H293T_high_vsn <- colMeans(model_vsn[se_UVC$celltype == "H293T" & se_UVC$section == "high", ])
H293T_med_vsn <- colMeans(model_vsn[se_UVC$celltype == "H293T" & se_UVC$section == "med", ])
H293T_low_vsn <- colMeans(model_vsn[se_UVC$celltype == "H293T" & se_UVC$section == "low", ])
HepG2_high_vsn <- colMeans(model_vsn[se_UVC$celltype == "HepG2" & se_UVC$section == "high", ])
HepG2_med_vsn <- colMeans(model_vsn[se_UVC$celltype == "HepG2" & se_UVC$section == "med", ])
HepG2_low_vsn <- colMeans(model_vsn[se_UVC$celltype == "HepG2" & se_UVC$section == "low", ])

#Fit the data
fit1_norm_vsn = lmFit(assay(se_UVC), design = model.matrix(~ celltype * section, colData(se_UVC)))

#Get contrasts
contrast_matrix_vsn <- cbind("H293T_high.H293T_low" = H293T_high_vsn - H293T_low_vsn, "H293T_med.H293T_low" = H293T_med_vsn - H293T_low_vsn,
                             "HepG2_high.HepG2_low" = HepG2_high_vsn - HepG2_low_vsn, "HepG2_med.HepG2_low" = HepG2_med_vsn - HepG2_low_vsn,
                             "HepG2_high.H293T_high" = HepG2_high_vsn - H293T_high_vsn, "HepG2_med.H293T_med" = HepG2_med_vsn - H293T_med_vsn,
                             "HepG2_low.H293T_low" = HepG2_low_vsn - H293T_low_vsn)

#Second fit
fit2_norm_vsn <- contrasts.fit(fit1_norm_vsn, contrasts = contrast_matrix_vsn)
fit3_norm_vsn <- eBayes(fit2_norm_vsn)

#Get results
res_norm_vsn <- extract_results_table(fit3_norm_vsn) 

#Interaction analysis
fit1_norm_int_vsn = lmFit(assay(se_UVC), design = model.matrix(~ celltype * section, colData(se_UVC)))
fit2_norm_int_vsn <- eBayes(fit1_norm_int_vsn)
int_norm_vsn <- topTable(fit2_norm_int_vsn, coef = c("celltypeHepG2:sectionlow", "celltypeHepG2:sectionmed"), number = length(rownames(se_UVC))) 
int_norm_vsn$gene <- rownames(int_norm_vsn)

#Get significance labels
res_norm_vsn_HL <- merge(res_norm_vsn$H293T_high.H293T_low, res_norm_vsn$HepG2_high.HepG2_low, by="gene")
res_norm_vsn_HL <- merge( res_norm_vsn_HL,int_norm_vsn,  by = "gene")
res_norm_vsn_HL$int_sign <- res_norm_vsn_HL$adj.P.Val < 0.1

res_norm_vsn_ML <- merge(res_norm_vsn$H293T_med.H293T_low, res_norm_vsn$HepG2_med.HepG2_low, by="gene")
res_norm_vsn_ML <- merge( res_norm_vsn_ML,int_norm_vsn,  by = "gene")
res_norm_vsn_ML$int_sign <- res_norm_vsn_ML$adj.P.Val < 0.1

res_norm_vsn_HL$order_sign <- paste(res_norm_vsn_HL$logFC.x > 0 & res_norm_vsn_HL$fdr.x < 0.05, #293T up
                                           res_norm_vsn_HL$logFC.y > 0 & res_norm_vsn_HL$fdr.y < 0.05, #HepG2 up
                                           res_norm_vsn_HL$logFC.x < 0 & res_norm_vsn_HL$fdr.x < 0.05, #293T down
                                           res_norm_vsn_HL$logFC.y < 0 & res_norm_vsn_HL$fdr.y < 0.05, #HepG2 down
                                           sep = "_")

res_norm_vsn_ML$order_sign <- paste(res_norm_vsn_ML$logFC.x > 0 & res_norm_vsn_ML$fdr.x < 0.05, #293T up
                                           res_norm_vsn_ML$logFC.y > 0 & res_norm_vsn_ML$fdr.y < 0.05, #HepG2 up
                                           res_norm_vsn_ML$logFC.x < 0 & res_norm_vsn_ML$fdr.x < 0.05, #293T down
                                           res_norm_vsn_ML$logFC.y < 0 & res_norm_vsn_ML$fdr.y < 0.05, #HepG2 down
                                           sep = "_")
resHL <- res_norm_vsn_HL
colnames(resHL)[-1] <- paste("HL", colnames(resHL)[-1], sep="_")
resML <- res_norm_vsn_ML
colnames(resML)[-1] <- paste("ML", colnames(resML)[-1], sep="_")
res_evr <- merge(resHL, resML, by = "gene", all.x = TRUE, all.y = TRUE)

#Keep only UVC proteins that are in both Label-free and TMT experiments
bulk_HEK293T <- read.delim(paste("/Users/lducoli/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/Bulk_analysis/2_UVC_enriched/", rbp , "_HEK293T_clipper_results.txt", sep = ""), header = TRUE, row.names = 1)
bulk_HEK293T <- subset(bulk_HEK293T, FDR < 0.1 & logFC > log2(3))

bulk_HepG2 <- read.delim(paste("/Users/lducoli/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/Bulk_analysis/2_UVC_enriched/", rbp , "_HepG2_clipper_results.txt", sep = ""), header = TRUE, row.names = 1)
bulk_HepG2 <- subset(bulk_HepG2, FDR < 0.1 & logFC > log2(3))

bulk_UVC_gene <- unique(c(rownames(bulk_HepG2), rownames(bulk_HEK293T)))
TMT_UVC_gene <- get(paste(rbp2, "_DEP_UVC", sep = ""))$all$name[get(paste(rbp2, "_DEP_UVC", sep = ""))$all$UVC == "UVC"]

intersect <- intersect(bulk_UVC_gene, TMT_UVC_gene)

res_norm_vsn_HL <- res_norm_vsn_HL[res_norm_vsn_HL$gene %in%  intersect,]
res_norm_vsn_ML <- res_norm_vsn_ML[res_norm_vsn_ML$gene %in%  intersect,]

#Scatter plot
vsn_max_axis <- c(max(c( max(res_norm_vsn_HL$logFC.x, res_norm_vsn_HL$logFC.y),max(res_norm_vsn_ML$logFC.x, res_norm_vsn_ML$logFC.y))))
vsn_min_axis <- c(max(abs(min(res_norm_vsn_HL$logFC.x, res_norm_vsn_HL$logFC.y)),abs(min(res_norm_vsn_ML$logFC.x, res_norm_vsn_ML$logFC.y))))

colors <- c("FALSE" = "#868686","TRUE" =  "#ff8d2a")

colors2 <- c("FALSE_FALSE_FALSE_FALSE" = "#868686", "FALSE_TRUE_TRUE_FALSE" = "#770003", "TRUE_FALSE_FALSE_TRUE" = "#770003", "FALSE_FALSE_FALSE_TRUE" = "#2256ca" ,"FALSE_FALSE_TRUE_FALSE" = "#2256ca",
                    "FALSE_FALSE_TRUE_TRUE" = "#2256ca", "FALSE_TRUE_FALSE_FALSE" = "#00a04c", "TRUE_FALSE_FALSE_FALSE" = "#00a04c",
                    "TRUE_TRUE_FALSE_FALSE" =  "#00a04c")

res_norm_vsn_HL$gel <- "High"
res_norm_vsn_ML$gel <- "Medium"

res_norm_vsn_all <- rbind(res_norm_vsn_HL, res_norm_vsn_ML)

all.ggplot <- ggplot(data=res_norm_vsn_all, aes(x=logFC.x, y=logFC.y)) + geom_vline(xintercept = 0) + geom_hline(yintercept = 0) +  geom_abline(intercept = 0, linetype=2) + 
  geom_point(shape=21, size=2, aes(col=int_sign, fill = order_sign)) + #aes(col = int_sign)) +
  labs(title = paste(rbp, " (", nrow(res_norm_vsn_all)/2, " RDAPs)" , sep = "")  , x = expression("log2FC high vs. low sections in H293T"), y =
         expression("log2FC high vs. low sections in HepG2")) +
  scale_color_manual(values = colors) +
  scale_fill_manual(values = colors2) +
  ggrepel::geom_text_repel(data = res_norm_vsn_all[res_norm_vsn_all$int_sign == "TRUE",], aes(label = gene), size = 3, box.padding = unit(0.1, "lines"),
                           point.padding = unit(0.1, "lines"), segment.size = 0.5,max.overlaps = Inf) + 
  theme_bw() +
  theme( panel.grid.major = element_blank(), legend.position = "none",
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_blank(),
        plot.title = element_text(hjust = 0.5)) + facet_grid(~ gel) + xlim(-2.5,2.5)+ ylim(-2.5,2.5)

return(list(se_UVC = se_UVC, res_all = res_norm_vsn_all, res_evr = res_evr, plot_all = all.ggplot, UVC_genes = intersect))
}

ABCF1_DEP <- DEP_analysis(ABCF1_DEP_UVC$se, "ABCF1", "ABCF1" )
DDX5_DEP <- DEP_analysis(DDX5_DEP_UVC$se, "DDX5", "DDX5")
FUS_DEP <- DEP_analysis(FUS_DEP_UVC$se, "FUS", "FUS")
hnA2B1_DEP <- DEP_analysis(hnA2B1_DEP_UVC$se, "HNRNPA2B1", "hnA2B1")
hnC_DEP <- DEP_analysis(hnC_DEP_UVC$se, "HNRNPC", "hnC")
hnM_DEP <- DEP_analysis(hnM_DEP_UVC$se, "HNRNPM", "hnM")
hnU_DEP <- DEP_analysis(hnU_DEP_UVC$se, "HNRNPU", "hnU")
ILF2_DEP <- DEP_analysis(ILF2_DEP_UVC$se, "ILF2", "ILF2")
ILF3_DEP <- DEP_analysis(ILF3_DEP_UVC$se, "ILF3", "ILF3")
NAT10_DEP <- DEP_analysis(NAT10_DEP_UVC$se, "NAT10", "NAT10")
NONO_DEP <- DEP_analysis(NONO_DEP_UVC$se, "NONO", "NONO")
RBFOX2_DEP <- DEP_analysis(RBFOX2_DEP_UVC$se, "RBFOX2", "RBFOX2")
SFPQ_DEP <- DEP_analysis(SFPQ_DEP_UVC$se, "SFPQ", "SFPQ")

head(ABCF1_DEP$res_all, n = 2)
write.table(ABCF1_DEP$res_evr, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/ABCF1_TMT_twofactor_results.txt", row.names = FALSE, sep = "\t")
write.table(DDX5_DEP$res_evr, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/DDX5_TMT_twofactor_results.txt", row.names = FALSE, sep = "\t")
write.table(FUS_DEP$res_evr, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/FUS_TMT_twofactor_results.txt", row.names = FALSE, sep = "\t")
write.table(hnA2B1_DEP$res_evr, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/hnA2B1_TMT_twofactor_results.txt", row.names = FALSE, sep = "\t")
write.table(hnC_DEP$res_evr, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/hnC_TMT_twofactor_results.txt", row.names = FALSE, sep = "\t")
write.table(hnM_DEP$res_evr, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/hnM_TMT_twofactor_results.txt", row.names = FALSE, sep = "\t")
write.table(hnU_DEP$res_evr, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/hnU_TMT_twofactor_results.txt", row.names = FALSE, sep = "\t")
write.table(ILF2_DEP$res_evr, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/ILF2_TMT_twofactor_results.txt", row.names = FALSE, sep = "\t")
write.table(ILF3_DEP$res_evr, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/ILF3_TMT_twofactor_results.txt", row.names = FALSE, sep = "\t")
write.table(NAT10_DEP$res_evr, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/NAT10_TMT_twofactor_results.txt", row.names = FALSE, sep = "\t")
write.table(NONO_DEP$res_evr, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/NONO_TMT_twofactor_results.txt", row.names = FALSE, sep = "\t")
write.table(RBFOX2_DEP$res_evr, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/RBFOX2_TMT_twofactor_results.txt", row.names = FALSE, sep = "\t")
write.table(SFPQ_DEP$res_evr, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/SFPQ_TMT_twofactor_results.txt", row.names = FALSE, sep = "\t")

write.table(ABCF1_DEP$res_all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/ABCF1_TMT_twofactor_results_sub.txt", row.names = FALSE, sep = "\t")
write.table(DDX5_DEP$res_all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/DDX5_TMT_twofactor_results_sub.txt", row.names = FALSE, sep = "\t")
write.table(FUS_DEP$res_all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/FUS_TMT_twofactor_results_sub.txt", row.names = FALSE, sep = "\t")
write.table(hnA2B1_DEP$res_all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/hnA2B1_TMT_twofactor_results_sub.txt", row.names = FALSE, sep = "\t")
write.table(hnC_DEP$res_all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/hnC_TMT_twofactor_results_sub.txt", row.names = FALSE, sep = "\t")
write.table(hnM_DEP$res_all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/hnM_TMT_twofactor_results_sub.txt", row.names = FALSE, sep = "\t")
write.table(hnU_DEP$res_all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/hnU_TMT_twofactor_results_sub.txt", row.names = FALSE, sep = "\t")
write.table(ILF2_DEP$res_all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/ILF2_TMT_twofactor_results_sub.txt", row.names = FALSE, sep = "\t")
write.table(ILF3_DEP$res_all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/ILF3_TMT_twofactor_results_sub.txt", row.names = FALSE, sep = "\t")
write.table(NAT10_DEP$res_all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/NAT10_TMT_twofactor_results_sub.txt", row.names = FALSE, sep = "\t")
write.table(NONO_DEP$res_all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/NONO_TMT_twofactor_results_sub.txt", row.names = FALSE, sep = "\t")
write.table(RBFOX2_DEP$res_all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/RBFOX2_TMT_twofactor_results_sub.txt", row.names = FALSE, sep = "\t")
write.table(SFPQ_DEP$res_all, file = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/3_DEP/SFPQ_TMT_twofactor_results_sub.txt", row.names = FALSE, sep = "\t")

4. Visualization

Scatter plot of logFC high vs. low

We plot the results of the two factor analysis. We displayed only proteins that were categorized as RDAPs (UVC-enriched) in both label-free and TMT datasets. Green dots: High-MW-zone RDAPs, Blue dots: Low-MW-zone RDAPs, red dots: ambivalent RDAPs, orange border: FDR < 0.1 between HEK293T and HepG2.

ggarrange(ABCF1_DEP$plot_all,DDX5_DEP$plot_all,FUS_DEP$plot_all,hnA2B1_DEP$plot_all,
          hnC_DEP$plot_all,hnM_DEP$plot_all,hnU_DEP$plot_all,ILF2_DEP$plot_all,
          ILF3_DEP$plot_all,NAT10_DEP$plot_all,NONO_DEP$plot_all,RBFOX2_DEP$plot_all,
          SFPQ_DEP$plot_all, ncol = 3, nrow = 5)

# Save the bubble plot as pdf
pdf("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/4_Visualization/DEP_results_twofactors.pdf", height = 20, width = 22)
ggarrange(ABCF1_DEP$plot_all,DDX5_DEP$plot_all,FUS_DEP$plot_all,hnA2B1_DEP$plot_all,
          hnC_DEP$plot_all,hnM_DEP$plot_all,hnU_DEP$plot_all,ILF2_DEP$plot_all,
          ILF3_DEP$plot_all,NAT10_DEP$plot_all,NONO_DEP$plot_all,RBFOX2_DEP$plot_all,
          SFPQ_DEP$plot_all, ncol = 3, nrow = 5)
dev.off()

CDF of high-MW-zone and low-MW-zone

#Load MW table
MW_prot <- read.delim("/Users/lducoli/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/0_TMT_data/Tot_prot_MW.txt", header = TRUE)
MW_prot$gene <- gsub("\\..*","",MW_prot$gene)
MW_prot <- unique(MW_prot)

#Group definition
H293T_lvls <- list(none = "FALSE_FALSE", Low_MW_zone = "FALSE_TRUE", High_MW_zone = "TRUE_FALSE")
HepG2_lvls <- list(none = "FALSE_FALSE", Low_MW_zone = "FALSE_TRUE", High_MW_zone = "TRUE_FALSE")

#Determine the class for each cell type
get_MW <- function(rbp) {
  data <- get(paste(rbp, "_DEP", sep = ""))
  data$res_all$HEK293T_order_sign <- factor(paste(data$res_all$logFC.x > 0 & data$res_all$fdr.x < 0.05, 
                                                  data$res_all$logFC.x < 0 & data$res_all$fdr.x < 0.05, sep = "_"))
  levels(data$res_all$HEK293T_order_sign) <- H293T_lvls
  data$res_all$HepG2_order_sign <- factor(paste(data$res_all$logFC.y > 0 & data$res_all$fdr.y < 0.05, 
                                                data$res_all$logFC.y < 0 & data$res_all$fdr.y < 0.05, sep = "_"))
  levels(data$res_all$HepG2_order_sign) <- HepG2_lvls
  colnames(data$res_all)[1] <- "name"
  MW <- merge(data$res_all, get(rbp)[,1:2], by = c("name"), all.x = TRUE)
  MW <- merge(MW, MW_prot, by = c("ID"))
  return(list(data = data$res_all, MW = MW))
}

ABCF1_MW <- get_MW("ABCF1" )
DDX5_MW <- get_MW("DDX5")
FUS_MW <- get_MW("FUS")
hnA2B1_MW <- get_MW("hnA2B1")
hnC_MW <- get_MW("hnC")
hnM_MW <- get_MW("hnM")
hnU_MW <- get_MW("hnU")
ILF2_MW <- get_MW("ILF2")
ILF3_MW <- get_MW("ILF3")
NAT10_MW <- get_MW("NAT10")
NONO_MW <- get_MW("NONO")
RBFOX2_MW <- get_MW("RBFOX2")
SFPQ_MW <- get_MW("SFPQ")

#Get the proteins
Tot_UVC <- unique(rbind(unique(ABCF1_MW$MW[,c(1,2,20,21,22)]), unique(DDX5_MW$MW[,c(1,2,20,21,22)]), unique(FUS_MW$MW[,c(1,2,20,21,22)]),
                        unique(hnA2B1_MW$MW[,c(1,2,20,21,22)]), unique(hnC_MW$MW[,c(1,2,20,21,22)]), unique(hnM_MW$MW[,c(1,2,20,21,22)]),
                        unique(hnU_MW$MW[,c(1,2,20,21,22)]), unique(ILF2_MW$MW[,c(1,2,20,21,22)]), unique(ILF3_MW$MW[,c(1,2,20,21,22)]),
                        unique(NAT10_MW$MW[,c(1,2,20,21,22)]),unique(NONO_MW$MW[,c(1,2,20,21,22)]),unique(RBFOX2_MW$MW[,c(1,2,20,21,22)]), 
                        unique(SFPQ_MW$MW[,c(1,2,20,21,22)])))


Tot_UVC_HEK293T <- Tot_UVC[,c(1:3,5)]
Tot_UVC_HEK293T$celltype <- "HEK293T"
colnames(Tot_UVC_HEK293T)[3] <- "order_sign"
Tot_UVC_HEK293T <- unique(Tot_UVC_HEK293T)

Tot_UVC_HepG2 <- Tot_UVC[,c(1:2,4,5)]
Tot_UVC_HepG2$celltype <- "HepG2"
colnames(Tot_UVC_HepG2)[3] <- "order_sign"
Tot_UVC_HepG2 <- unique(Tot_UVC_HepG2)

Tot_UVC <- rbind(Tot_UVC_HEK293T, Tot_UVC_HepG2)

ggplot(data=Tot_UVC, aes(x=MW, group=order_sign, colour=order_sign)) +
  stat_ecdf() +
  ggtitle("Cumulative fraction of high and low-MW RDAP molecular weights (MW)") +
  xlim(0, 400) +
  scale_color_manual(values = c("#1d1d1b", "#2256ca", "#00a04c")) +
  geom_vline(xintercept = c(60,350), linetype="dashed", color = "black", size=0.5) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_blank(),
        plot.title = element_text(hjust = 0.5)) + facet_grid(~celltype)

# Save the bubble plot as pdf
pdf("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/4_Visualization/CDF_high_low_MW.pdf", height = 5, width = 10)
ggplot(data=Tot_UVC, aes(x=MW, group=order_sign, colour=order_sign)) +
  stat_ecdf() +
  ggtitle("Cumulative fraction of high and low-MW RDAP molecular weights (MW)") +
  xlim(0, 400) +
  scale_color_manual(values = c("#1d1d1b",  "#2256ca", "#00a04c")) +
  geom_vline(xintercept = c(60,350), linetype="dashed", color = "black", size=0.5) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_blank(),
        plot.title = element_text(hjust = 0.5)) + facet_grid(~celltype)
dev.off()
#KS test
data.frame( row.names = c("HEK293T", "HepG2"), 
            highvsbackground = c(ks.test(Tot_UVC_HEK293T$MW[Tot_UVC_HEK293T$order_sign == "none"], 
                                         Tot_UVC_HEK293T$MW[Tot_UVC_HEK293T$order_sign == "High_MW_zone"] , alternative = 'two.sided')$p.value,
                                 ks.test(Tot_UVC_HepG2$MW[Tot_UVC_HepG2$order_sign == "none"], 
                                         Tot_UVC_HepG2$MW[Tot_UVC_HepG2$order_sign == "High_MW_zone"] , alternative = 'two.sided')$p.value),
            lowvsbackground = c(ks.test(Tot_UVC_HEK293T$MW[Tot_UVC_HEK293T$order_sign == "none"], 
                                         Tot_UVC_HEK293T$MW[Tot_UVC_HEK293T$order_sign == "Low_MW_zone"] , alternative = 'two.sided')$p.value,
                                 ks.test(Tot_UVC_HepG2$MW[Tot_UVC_HepG2$order_sign == "none"], 
                                         Tot_UVC_HepG2$MW[Tot_UVC_HepG2$order_sign == "Low_MW_zone"] , alternative = 'two.sided')$p.value),
            highvslow = c(ks.test(Tot_UVC_HEK293T$MW[Tot_UVC_HEK293T$order_sign == "Low_MW_zone"], 
                                         Tot_UVC_HEK293T$MW[Tot_UVC_HEK293T$order_sign == "High_MW_zone"] , alternative = 'two.sided')$p.value,
                                 ks.test(Tot_UVC_HepG2$MW[Tot_UVC_HepG2$order_sign == "Low_MW_zone"], 
                                         Tot_UVC_HepG2$MW[Tot_UVC_HepG2$order_sign == "High_MW_zone"] , alternative = 'two.sided')$p.value))

Heatmap of slope values

We calculated the slope from normalized log2-transformed TMT intensities across the three RNP subzones. We then displayed these values as an heatmap and through clustering, we have identified 20 RDAPs with not only multi-protein assembly tendency but also cell-type differences.

merge.all <- function(x, ..., by = "row.names") {
  L <- list(...)
  for (i in seq_along(L)) {
    x <- merge(x, L[[i]], by = by, all = TRUE)
    rownames(x) <- x$Row.names
    x$Row.names <- NULL
  }
  return(x)
}

#Get all proteins with cell-type enrichment
Tot_TMT_int <- unique(c(unique(ABCF1_DEP$res_all$gene[ABCF1_DEP$res_all$int_sign == "TRUE"]),
                        unique(DDX5_DEP$res_all$gene[DDX5_DEP$res_all$int_sign == "TRUE"]),
                        unique(FUS_DEP$res_all$gene[FUS_DEP$res_all$int_sign == "TRUE"]),
                        unique(hnA2B1_DEP$res_all$gene[hnA2B1_DEP$res_all$int_sign == "TRUE"]),
                        unique(hnC_DEP$res_all$gene[hnC_DEP$res_all$int_sign == "TRUE"]),
                        unique(hnM_DEP$res_all$gene[hnM_DEP$res_all$int_sign == "TRUE"]),
                        unique(hnU_DEP$res_all$gene[hnU_DEP$res_all$int_sign == "TRUE"]),
                        unique(ILF2_DEP$res_all$gene[ILF2_DEP$res_all$int_sign == "TRUE"]),
                        unique(ILF3_DEP$res_all$gene[ILF3_DEP$res_all$int_sign == "TRUE"]),
                        unique(NAT10_DEP$res_all$gene[NAT10_DEP$res_all$int_sign == "TRUE"]),
                        unique(NONO_DEP$res_all$gene[NONO_DEP$res_all$int_sign == "TRUE"]),
                        unique(RBFOX2_DEP$res_all$gene[RBFOX2_DEP$res_all$int_sign == "TRUE"]),
                        unique(SFPQ_DEP$res_all$gene[SFPQ_DEP$res_all$int_sign == "TRUE"])))

#Average the replicates
rowMeans_function <- function(data, rbp) {
data <- as.data.frame(assay(data))
data$H293T_low_avg <- rowMeans(data[,c(1,4)])
data$H293T_med_avg <- rowMeans(data[,c(2,5)])
data$H293T_high_avg <- rowMeans(data[,c(3,6)])
data$HepG2_low_avg <- rowMeans(data[,c(7,10)])
data$HepG2_med_avg <- rowMeans(data[,c(8,11)])
data$HepG2_high_avg <- rowMeans(data[,c(9,12)])
data <- data[,13:18]
colnames(data) <- paste(rbp, colnames(data), sep="_")
data2 <- data[rownames(data) %in% Tot_TMT_int,]
data2 <- data2[rownames(data2) %in% get(paste(rbp, "_DEP", sep = ""))$UVC_genes,]
return(list(unflt = data, flt = data2))
}

ABCF1_data <- rowMeans_function(ABCF1_DEP$se_UVC, "ABCF1")
DDX5_data <- rowMeans_function(DDX5_DEP$se_UVC, "DDX5")
FUS_data <- rowMeans_function(FUS_DEP$se_UVC, "FUS")
hnA2B1_data <- rowMeans_function(hnA2B1_DEP$se_UVC, "hnA2B1")
hnC_data <- rowMeans_function(hnC_DEP$se_UVC, "hnC")
hnM_data <- rowMeans_function(hnM_DEP$se_UVC, "hnM")
hnU_data <- rowMeans_function(hnU_DEP$se_UVC, "hnU")
ILF2_data <- rowMeans_function(ILF2_DEP$se_UVC, "ILF2")
ILF3_data <- rowMeans_function(ILF3_DEP$se_UVC, "ILF3")
NAT10_data <- rowMeans_function(NAT10_DEP$se_UVC, "NAT10")
NONO_data <- rowMeans_function(NONO_DEP$se_UVC, "NONO")
RBFOX2_data <- rowMeans_function(RBFOX2_DEP$se_UVC, "RBFOX2")
SFPQ_data <- rowMeans_function(SFPQ_DEP$se_UVC, "SFPQ")


all_data <- merge.all(ABCF1_data$flt,DDX5_data$flt,FUS_data$flt,hnA2B1_data$flt,hnC_data$flt,hnM_data$flt,hnU_data$flt,ILF2_data$flt,ILF3_data$flt,NAT10_data$flt,NONO_data$flt,RBFOX2_data$flt,SFPQ_data$flt)

l <- list(c(grep("ABCF1_H293T", colnames(all_data)),grep("ABCF1_HepG2", colnames(all_data))),
          c(grep("DDX5_H293T", colnames(all_data)),grep("DDX5_HepG2", colnames(all_data))),
          c(grep("FUS_H293T", colnames(all_data)),grep("FUS_HepG2", colnames(all_data))),
          c(grep("hnA2B1_H293T", colnames(all_data)),grep("hnA2B1_HepG2", colnames(all_data))),
          c(grep("hnC_H293T", colnames(all_data)),grep("hnC_HepG2", colnames(all_data))),
          c(grep("hnM_H293T", colnames(all_data)),grep("hnM_HepG2", colnames(all_data))),
          c(grep("hnU_H293T", colnames(all_data)),grep("hnU_HepG2", colnames(all_data))),
          c(grep("ILF2_H293T", colnames(all_data)),grep("ILF2_HepG2", colnames(all_data))),
          c(grep("ILF3_H293T", colnames(all_data)),grep("ILF3_HepG2", colnames(all_data))),
          c(grep("NAT10_H293T", colnames(all_data)),grep("NAT10_HepG2", colnames(all_data))),
          c(grep("NONO_H293T", colnames(all_data)),grep("NONO_HepG2", colnames(all_data))),
          c(grep("RBFOX2_H293T", colnames(all_data)),grep("RBFOX2_HepG2", colnames(all_data))),
          c(grep("SFPQ_H293T", colnames(all_data)),grep("SFPQ_HepG2", colnames(all_data)))
          )

for (i in 1:length(l)) {
  all_data[,l[[i]]] <- t(scale(t(all_data[,l[[i]]]), scale = FALSE, center = TRUE))
}

slope  <-  function(x){
  if(all(is.na(x)))
    return(NA)
  else
    return(coef(lm(x~I(1:3)))[2])
}

#Slope
all_data$ABCF1_H293T_slope <- apply(all_data[,grep("ABCF1_H293T", colnames(all_data))], 1,slope)
all_data$ABCF1_HepG2_slope <- apply(all_data[,grep("ABCF1_HepG2", colnames(all_data))], 1,slope)
all_data$DDX5_H293T_slope <- apply(all_data[,grep("DDX5_H293T", colnames(all_data))], 1,slope)
all_data$DDX5_HepG2_slope <- apply(all_data[,grep("DDX5_HepG2", colnames(all_data))], 1,slope)
all_data$FUS_H293T_slope <- apply(all_data[,grep("FUS_H293T", colnames(all_data))], 1,slope)
all_data$FUS_HepG2_slope <- apply(all_data[,grep("FUS_HepG2", colnames(all_data))], 1,slope)
all_data$hnA2B1_H293T_slope <- apply(all_data[,grep("hnA2B1_H293T", colnames(all_data))], 1,slope)
all_data$hnA2B1_HepG2_slope <- apply(all_data[,grep("hnA2B1_HepG2", colnames(all_data))], 1,slope)
all_data$hnC_H293T_slope <- apply(all_data[,grep("hnC_H293T", colnames(all_data))], 1,slope)
all_data$hnC_HepG2_slope <- apply(all_data[,grep("hnC_HepG2", colnames(all_data))], 1,slope)
all_data$hnM_H293T_slope <- apply(all_data[,grep("hnM_H293T", colnames(all_data))], 1,slope)
all_data$hnM_HepG2_slope <- apply(all_data[,grep("hnM_HepG2", colnames(all_data))], 1,slope)
all_data$hnU_H293T_slope <- apply(all_data[,grep("hnU_H293T", colnames(all_data))], 1,slope)
all_data$hnU_HepG2_slope <- apply(all_data[,grep("hnU_HepG2", colnames(all_data))], 1,slope)
all_data$ILF2_H293T_slope <- apply(all_data[,grep("ILF2_H293T", colnames(all_data))], 1,slope)
all_data$ILF2_HepG2_slope <- apply(all_data[,grep("ILF2_HepG2", colnames(all_data))], 1,slope)
all_data$ILF3_H293T_slope <- apply(all_data[,grep("ILF3_H293T", colnames(all_data))], 1,slope)
all_data$ILF3_HepG2_slope <- apply(all_data[,grep("ILF3_HepG2", colnames(all_data))], 1,slope)
all_data$NAT10_H293T_slope <- apply(all_data[,grep("NAT10_H293T", colnames(all_data))], 1,slope)
all_data$NAT10_HepG2_slope <- apply(all_data[,grep("NAT10_HepG2", colnames(all_data))], 1,slope)
all_data$NONO_H293T_slope <- apply(all_data[,grep("NONO_H293T", colnames(all_data))], 1,slope)
all_data$NONO_HepG2_slope <- apply(all_data[,grep("NONO_HepG2", colnames(all_data))], 1,slope)
all_data$RBFOX2_H293T_slope <- apply(all_data[,grep("RBFOX2_H293T", colnames(all_data))], 1,slope)
all_data$RBFOX2_HepG2_slope <- apply(all_data[,grep("RBFOX2_HepG2", colnames(all_data))], 1,slope)
all_data$SFPQ_H293T_slope <- apply(all_data[,grep("SFPQ_H293T", colnames(all_data))], 1,slope)
all_data$SFPQ_HepG2_slope <- apply(all_data[,grep("SFPQ_HepG2", colnames(all_data))], 1,slope)

#Calculate the clustering
all_data_slope <- all_data[,grep("slope", colnames(all_data))]
all_data_slope_cl <- all_data_slope
all_data_slope_cl[is.na(all_data_slope_cl)] <- 0

hc <- hclust(dist(all_data_slope_cl), method = "ward.D2")

# Plot dendogram
plot(hclust(dist(all_data_slope_cl), method = "ward.D2"), hang = -1, cex=0.5)

# Save dendogram as pdf
pdf("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/4_Visualization/Dendo_slope.pdf")
plot(hclust(dist(all_data_slope_cl), method = "ward.D2"), hang = -1, cex=0.5)
dev.off()
#Annotation
annotation <- MW_prot[,1:2]
rownames(annotation) <- annotation$gene
annotation$gene <- NULL
annotation <- subset(annotation, rownames(annotation) %in% rownames(all_data))
annotation2 <- data.frame(row.names = colnames(all_data_slope), celltype = rep(c("HEK293T", "HepG2"), 13))
ann_colors <- list(MW = c(paletteer_c("grDevices::Peach", 30)), celltype =  c(HEK293T = "#f79516", HepG2 = "#00ace6"))

fromList <- function (input) {
  elements <- unique(unlist(input))
  data <- unlist(lapply(input, function(x) {
    x <- as.vector(match(elements, x))
  }))
  data[is.na(data)] <- as.integer(0)
  data[data != 0] <- as.integer(1)
  data <- data.frame(matrix(data, ncol = length(input), byrow = F))
  data <- data[which(rowSums(data) != 0), ]
  names(data) <- names(input)
  row.names(data) <- elements
  return(data)
}

#Heatmap color range
my.breaks <- c(seq(-0.75, 0.75, by=0.01))
my.colors <- rev(c(paletteer_c("ggthemes::Green-Blue-White Diverging", length(my.breaks))))

lt.tsk = list(ABCF1 = unique(ABCF1_DEP$res_all$gene[ABCF1_DEP$res_all$int_sign == "TRUE"]),
              DDX5 = unique(DDX5_DEP$res_all$gene[DDX5_DEP$res_all$int_sign == "TRUE"]),
              FUS = unique(FUS_DEP$res_all$gene[FUS_DEP$res_all$int_sign == "TRUE"]),
              hnA2B1 = unique(hnA2B1_DEP$res_all$gene[hnA2B1_DEP$res_all$int_sign == "TRUE"]),
              hnC = unique(hnC_DEP$res_all$gene[hnC_DEP$res_all$int_sign == "TRUE"]),
              hnM = unique(hnM_DEP$res_all$gene[hnM_DEP$res_all$int_sign == "TRUE"]),
              hnU = unique(hnU_DEP$res_all$gene[hnU_DEP$res_all$int_sign == "TRUE"]),
              ILF2 = unique(ILF2_DEP$res_all$gene[ILF2_DEP$res_all$int_sign == "TRUE"]),
              ILF3 = unique(ILF3_DEP$res_all$gene[ILF3_DEP$res_all$int_sign == "TRUE"]),
              NAT10 = unique(NAT10_DEP$res_all$gene[NAT10_DEP$res_all$int_sign == "TRUE"]),
              NONO = unique(NONO_DEP$res_all$gene[NONO_DEP$res_all$int_sign == "TRUE"]),
              RBFOX2 = unique(RBFOX2_DEP$res_all$gene[RBFOX2_DEP$res_all$int_sign == "TRUE"]),
              SFPQ = unique(SFPQ_DEP$res_all$gene[SFPQ_DEP$res_all$int_sign == "TRUE"])
              )

# Binary table with colnames:
sign.proteins <- fromList(lt.tsk)

labels <- sign.proteins[,c(rep(1:13, each = 2))]
colnames(labels) <- colnames(all_data[,grep("slope", colnames(all_data))])
labels[labels == 1] <- "*"
labels[labels == 0] <- ""
labels <- labels[match(rownames(all_data_slope), rownames(labels)),]

#Heatmap
pheatmap <- pheatmap(
  mat               = t(all_data_slope),
  annotation_col = annotation,
  annotation_row = annotation2,
  annotation_colors = ann_colors,
  cellwidth = 8,
  cellheight = 6,
  display_numbers = t(labels), 
  fontsize_number=5.5,
  color = my.colors,
  breaks = my.breaks,
  show_colnames     = TRUE,
  show_rownames     = TRUE,
  drop_levels       = TRUE,
  fontsize          = 5.5,
  cluster_rows      = FALSE,
  cluster_cols      = hc,
  gaps_row = c(2,4,6,8,10,12,14,16,18,20,22,24),
  na_col = "lightgrey"
)

# Save the heatmap
pheatmap(
  mat               = t(all_data_slope),
  annotation_col = annotation,
  annotation_row = annotation2,
  annotation_colors = ann_colors,
  cellwidth = 8,
  cellheight = 6,
  display_numbers = t(labels), 
  fontsize_number=5.5,
  color = my.colors,
  breaks = my.breaks,
  show_colnames     = TRUE,
  show_rownames     = TRUE,
  drop_levels       = TRUE,
  fontsize          = 5.5,
  cluster_rows      = FALSE,
  cluster_cols      = hc,
  gaps_row = c(2,4,6,8,10,12,14,16,18,20,22,24),
  na_col = "lightgrey",
  width = 10,
  height = 5,
  filename = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/4_Visualization/HM_TMT_sign_int.pdf"
)

Network analysis between HEK293T and HepG2

Network analysis was performed between the 13 RBPs and the high-MW RDAP cluster member identified above. We used a spin-glass model to identify communities in the different networks.

weight.community <- function(row,membership,weigth.within,weight.between){
  if(as.numeric(membership[which(names(membership)==row[1])])==as.numeric(membership[which(names(membership)==row[2])])){weight=weigth.within
  }else{
    weight=weight.between
    }
  return(weight)
  }

Get_network_table <- function(rbp) {
  gel_order <- get(paste(rbp, "_MW", sep = ""))$MW[,c(1,2,17,19,20,21,23)]
  gel_order <- subset(gel_order, name %in% cluster$name[cluster$cluster == 2])
  gel_order$rbp <- rbp
  return(gel_order)
}

#Cluster
cluster <- data.frame(name = rownames(all_data_slope), cluster = cutree(hc, k=3))
rownames(cluster) <- NULL

#get table for network
ABCF1_g <- Get_network_table("ABCF1")
DDX5_g <- Get_network_table("DDX5")
FUS_g <- Get_network_table("FUS")
hnA2B1_g <- Get_network_table("hnA2B1")
hnA2B1_g$rbp <- "HNRNPA2B1"
hnC_g <- Get_network_table("hnC")
hnC_g$rbp <- "HNRNPC"
hnM_g <- Get_network_table("hnM")
hnM_g$rbp <- "HNRNPM"
hnU_g <- Get_network_table("hnU")
hnU_g$rbp <- "HNRNPU"
ILF2_g <- Get_network_table("ILF2")
ILF3_g <- Get_network_table("ILF3")
NAT10_g <- Get_network_table("NAT10")
NONO_g <- Get_network_table("NONO")
RBFOX2_g <- Get_network_table("RBFOX2")
SFPQ_g <- Get_network_table("SFPQ")

#Get connection for every RBP
gel_order <- rbind(ABCF1_g,DDX5_g,FUS_g,hnA2B1_g,hnC_g,hnM_g,hnU_g,ILF2_g,ILF3_g,NAT10_g,NONO_g,RBFOX2_g,SFPQ_g)

#Prepare dataset for clustering
gel_order$HEK293T_edges <- paste(gel_order$HEK293T_order_sign, gel_order$int_sign, sep="_")
gel_order$HEK293T_edges <- factor(gel_order$HEK293T_edges)
levels(gel_order$HEK293T_edges) <- c("#00a04c", "#770003","#2256ca", "orange", "#868686","#868686")
gel_order$HepG2_edges <- paste(gel_order$HepG2_order_sign, gel_order$int_sign, sep="_")
gel_order$HepG2_edges <- factor(gel_order$HepG2_edges)
levels(gel_order$HepG2_edges) <- c("#00a04c", "#770003", "orange", "#868686","#868686")

network_object <- function(data, celltype) {
  set.seed(5)
  nodes <- data.frame(name=unique(c(data$rbp,data$name)),
                    type=c(rep("Bait", 13), rep("RDAP", length(unique(c(data$rbp,data$name)))-13)),
                    size=c(rep(15, 13), rep(3, length(unique(c(data$rbp,data$name)))-13)))
  edges <- data.frame(from= data$rbp,
                        to=data$name,
                        celltype=data[,colnames(data) == paste(celltype, "_edges", sep = "")])
  g <- graph.data.frame(edges, directed=TRUE, vertices=nodes)
  g <- simplify(g, remove.multiple = F, remove.loops = T) 
  E(g)$color <- E(g)$celltype
  clp <- cluster_spinglass(as.undirected(g))
  mod <- modularity(clp)
  V(g)$community <- clp$membership
  E(g)$weight=apply(get.edgelist(g),1,weight.community,membership(clp),10,1)
  g$layout=layout.fruchterman.reingold(g,weights=E(g)$weight)
  colrs <- c(paletteer_dynamic("cartography::multi.pal", length(unique(clp$membership))))
  return(list(clp = clp, colrs = colrs, g = g, mod = mod))
}

#High
gel_order_high <- subset(gel_order, gel %in% "High")

gel_order_high_HEK293T <- subset(gel_order_high, HEK293T_edges != "#868686")
gel_order_high_HepG2 <- subset(gel_order_high, HepG2_edges != "#868686")

HEK293T_high <- network_object(gel_order_high_HEK293T, "HEK293T")
HepG2_high <- network_object(gel_order_high_HepG2, "HepG2")

#Medium
gel_order_medium <- subset(gel_order, gel %in% "Medium")

gel_order_medium_HEK293T <- subset(gel_order_medium, HEK293T_edges != "#868686")
gel_order_medium_HepG2 <- subset(gel_order_medium, HepG2_edges != "#868686")

HEK293T_medium <- network_object(gel_order_medium_HEK293T, "HEK293T")
HepG2_medium <- network_object(gel_order_medium_HepG2, "HepG2")

par(mfrow=c(1,4))
plot(HEK293T_high$g, edge.arrow.size=0.5, vertex.label.color="black", vertex.label.dist=0, vertex.color=HEK293T_high$colrs[V(HEK293T_high$g)$community], edge.curved=0.2)
plot(HepG2_high$g, edge.arrow.size=0.5, vertex.label.color="black", vertex.label.dist=0, vertex.color=HepG2_high$colrs[V(HepG2_high$g)$community], edge.curved=0.2)
plot(HEK293T_medium$g, edge.arrow.size=0.5, vertex.label.color="black", vertex.label.dist=0, vertex.color=HEK293T_medium$colrs[V(HEK293T_medium$g)$community], edge.curved=0.2)
plot(HepG2_medium$g, edge.arrow.size=0.5, vertex.label.color="black", vertex.label.dist=0, vertex.color=HepG2_medium$colrs[V(HepG2_medium$g)$community], edge.curved=0.2)

# Save dendogram as pdf
pdf("~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/4_Visualization/Network_highMW_RDAPs.pdf", width = 30, height = 10)
par(mfrow=c(1,4))
plot(HEK293T_high$g, edge.arrow.size=0.5, vertex.label.color="black", vertex.label.dist=0, vertex.color=HEK293T_high$colrs[V(HEK293T_high$g)$community], edge.curved=0.2)
plot(HepG2_high$g, edge.arrow.size=0.5, vertex.label.color="black", vertex.label.dist=0, vertex.color=HepG2_high$colrs[V(HepG2_high$g)$community], edge.curved=0.2)
plot(HEK293T_medium$g, edge.arrow.size=0.5, vertex.label.color="black", vertex.label.dist=0, vertex.color=HEK293T_medium$colrs[V(HEK293T_medium$g)$community], edge.curved=0.2)
plot(HepG2_medium$g, edge.arrow.size=0.5, vertex.label.color="black", vertex.label.dist=0, vertex.color=HepG2_medium$colrs[V(HepG2_medium$g)$community], edge.curved=0.2)
dev.off()

Reciprocal ratio Heatmaps for RBPs

Similar to what we have done with the label-free irCLIP-RNP data, we have generated a reciprocal heatmap that compare the ratio of intensities across the three RNP subzones for the 13 RBPs.

all_data2 <- merge.all(ABCF1_data$unflt,DDX5_data$unflt,FUS_data$unflt,hnA2B1_data$unflt,hnC_data$unflt,hnM_data$unflt,hnU_data$unflt,ILF2_data$unflt,ILF3_data$unflt,NAT10_data$unflt,NONO_data$unflt,RBFOX2_data$unflt,SFPQ_data$unflt)

get_ratio <- function(data, rbp) {
  data$H293T_rowsum <- rowSums(2^data[,grep(paste(rbp,"H293T", sep = "_"), colnames(data))])
  data$H293T_L_ratio <- 2^data[,grep(paste(rbp,"H293T_low_avg", sep = "_"), colnames(data))]/data$H293T_rowsum
  data$H293T_M_ratio <- 2^data[,grep(paste(rbp,"H293T_med_avg", sep = "_"), colnames(data))]/data$H293T_rowsum
  data$H293T_H_ratio <- 2^data[,grep(paste(rbp,"H293T_high_avg", sep = "_"), colnames(data))]/data$H293T_rowsum
  data$HepG2_rowsum <- rowSums(2^data[,grep(paste(rbp,"HepG2", sep = "_"), colnames(data))])
  data$HepG2_L_ratio <- 2^data[,grep(paste(rbp,"HepG2_low_avg", sep = "_"), colnames(data))]/data$HepG2_rowsum
  data$HepG2_M_ratio <- 2^data[,grep(paste(rbp,"HepG2_med_avg", sep = "_"), colnames(data))]/data$HepG2_rowsum
  data$HepG2_H_ratio <- 2^data[,grep(paste(rbp,"HepG2_high_avg", sep = "_"), colnames(data))]/data$HepG2_rowsum
  data$H293T_rowsum <- NULL
  data$HepG2_rowsum <- NULL
  colnames(data)[grep("ratio", colnames(data))] <- paste(rbp, colnames(data)[grep("ratio", colnames(data))], sep = "_")
  return(data)
}

ABCF1_ratio <- get_ratio(all_data2, "ABCF1")
DDX5_ratio <- get_ratio(all_data2, "DDX5")
FUS_ratio <- get_ratio(all_data2, "FUS")
hnA2B1_ratio <- get_ratio(all_data2, "hnA2B1")
hnC_ratio <- get_ratio(all_data2, "hnC")
hnM_ratio <- get_ratio(all_data2, "hnM")
hnU_ratio <- get_ratio(all_data2, "hnU")
ILF2_ratio <- get_ratio(all_data2, "ILF2")
ILF3_ratio <- get_ratio(all_data2, "ILF3")
NAT10_ratio <- get_ratio(all_data2, "NAT10")
NONO_ratio <- get_ratio(all_data2, "NONO")
RBFOX2_ratio <- get_ratio(all_data2, "RBFOX2")
SFPQ_ratio <- get_ratio(all_data2, "SFPQ")

rbp_list <- c("ABCF1", "DDX5", "FUS", "HNRNPA2B1", "HNRNPC", "HNRNPM", "HNRNPU", "ILF2", "ILF3", "NAT10", "NONO", "RBFOX2", "SFPQ")  

all_data_baitratio <- cbind(ABCF1_ratio[,grep("ratio", colnames(ABCF1_ratio))], 
                            DDX5_ratio[,grep("ratio", colnames(DDX5_ratio))],
                            FUS_ratio[,grep("ratio", colnames(FUS_ratio))],
                            hnA2B1_ratio[,grep("ratio", colnames(hnA2B1_ratio))],
                            hnC_ratio[,grep("ratio", colnames(hnC_ratio))],
                            hnM_ratio[,grep("ratio", colnames(hnM_ratio))],
                            hnU_ratio[,grep("ratio", colnames(hnU_ratio))],
                            ILF2_ratio[,grep("ratio", colnames(ILF2_ratio))],
                            ILF3_ratio[,grep("ratio", colnames(ILF3_ratio))],
                            NAT10_ratio[,grep("ratio", colnames(NAT10_ratio))],
                            NONO_ratio[,grep("ratio", colnames(NONO_ratio))],
                            RBFOX2_ratio[,grep("ratio", colnames(RBFOX2_ratio))],
                            SFPQ_ratio[,grep("ratio", colnames(SFPQ_ratio))])
all_data_baitratio <- all_data_baitratio[rbp_list,]

head(all_data_baitratio)
all_data_H293T <- all_data_baitratio[,grep("H293T", colnames(all_data_baitratio))]
colnames(all_data_H293T) <- gsub("_H293T", "", colnames(all_data_H293T))
rownames(all_data_H293T) <- paste(rownames(all_data_H293T), "H293T", sep = "_")
all_data_HepG2 <- all_data_baitratio[,grep("HepG2", colnames(all_data_baitratio))]
colnames(all_data_HepG2) <- gsub("_HepG2", "", colnames(all_data_HepG2))
rownames(all_data_HepG2) <- paste(rownames(all_data_HepG2), "HepG2", sep = "_")

all_data_baitratio2 <- rbind(all_data_H293T, all_data_HepG2)
#Pheatmap
my.breaks <- c(seq(0.2, 0.5, by=0.01)) 
my.colors <- c("#f9f9f9", rev(paletteer_c("grDevices::Greens 3" , length(my.breaks)-1)))

#Annotation
annotation_col <- data.frame(row.names = colnames(all_data_baitratio2), RNPzone = rep(c("RNPzone1", "RNPzone2", "RNPzone3"),13))
annotation_row <- data.frame(row.names = rownames(all_data_baitratio2), celltype = rep(c("HEK293T", "HepG2"),each = 13))

ann_colors2 <- list(RNPzone =  c(RNPzone1 = "#2256ca", RNPzone2 = "#b6de43", RNPzone3 = "#269541"), celltype =  c(HEK293T = "#f79516", HepG2 = "#00ace6"))

pheatmap <- pheatmap(
  mat               = all_data_baitratio2,
  annotation_col = annotation_col,
  annotation_row = annotation_row,
  annotation_colors = ann_colors2,
  cellwidth = 12,
  cellheight = 6,
  fontsize_number=5.5,
  color = my.colors,
  breaks = my.breaks,
  show_colnames     = TRUE,
  show_rownames     = TRUE,
  drop_levels       = TRUE,
  fontsize          = 5.5,
  cluster_rows      = FALSE,
  cluster_cols      = FALSE,
  gaps_col = c(3,6,9,12,15,18,21,24,27,30,33,36),
  gaps_row = c(13),
  na_col = "lightgrey"
)

# Save the heatmap
pheatmap(
  mat               = all_data_baitratio2,
  annotation_col = annotation_col,
  annotation_row = annotation_row,
  annotation_colors = ann_colors2,
  cellwidth = 12,
  cellheight = 6,
  fontsize_number=5.5,
  color = my.colors,
  breaks = my.breaks,
  show_colnames     = TRUE,
  show_rownames     = TRUE,
  drop_levels       = TRUE,
  fontsize          = 5.5,
  cluster_rows      = FALSE,
  cluster_cols      = FALSE,
  gaps_col = c(3,6,9,12,15,18,21,24,27,30,33,36),
  gaps_row = c(13),
  na_col = "lightgrey",
  width = 10,
  height = 3,
  filename = "~/Documents/Postdoc/PD_Projects/3_irCLIP-RNP/MS/TMT_analysis/4_Visualization/HM_TMT_reciprocal.pdf"
)

All the visualizations were saved as pdf and modified in illustrator.

sessionInfo()
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] rstatix_0.7.2               igraph_2.0.3               
##  [3] GGally_2.2.1                gprofiler2_0.2.3           
##  [5] ggpmisc_0.5.5               ggpp_0.5.6                 
##  [7] psych_2.4.3                 corrplot_0.92              
##  [9] paletteer_1.6.0             factoextra_1.0.7           
## [11] Clipper_0.0.0.9000          UpSetR_1.4.0               
## [13] ggpubr_0.6.0                DESeq2_1.38.3              
## [15] hexbin_1.28.3               viridis_0.6.5              
## [17] viridisLite_0.4.2           ggExtra_0.10.1             
## [19] textshape_1.7.3             pacman_0.5.1               
## [21] hrbrthemes_0.8.7            gplots_3.1.3.1             
## [23] RColorBrewer_1.1-3          pheatmap_1.0.12            
## [25] data.table_1.15.2           lubridate_1.9.3            
## [27] forcats_1.0.0               stringr_1.5.1              
## [29] dplyr_1.1.4                 purrr_1.0.2                
## [31] readr_2.1.5                 tidyr_1.3.1                
## [33] tibble_3.2.1                ggplot2_3.5.0              
## [35] tidyverse_2.0.0             DEP2_0.4.8.24              
## [37] R6_2.5.1                    limma_3.54.2               
## [39] MSnbase_2.24.2              ProtGenerics_1.30.0        
## [41] mzR_2.32.0                  Rcpp_1.0.12                
## [43] MsCoreUtils_1.10.0          SummarizedExperiment_1.28.0
## [45] Biobase_2.58.0              GenomicRanges_1.50.2       
## [47] GenomeInfoDb_1.34.9         IRanges_2.32.0             
## [49] S4Vectors_0.36.2            BiocGenerics_0.44.0        
## [51] MatrixGenerics_1.10.0       matrixStats_1.2.0          
## 
## loaded via a namespace (and not attached):
##   [1] SparseM_1.81                ggthemes_5.1.0             
##   [3] missForest_1.5              bit64_4.0.5                
##   [5] knitr_1.45                  DelayedArray_0.24.0        
##   [7] KEGGREST_1.38.0             RCurl_1.98-1.14            
##   [9] AnnotationFilter_1.22.0     doParallel_1.0.17          
##  [11] generics_0.1.3              preprocessCore_1.60.2      
##  [13] cowplot_1.1.3               RSQLite_2.3.5              
##  [15] proxy_0.4-27                bit_4.0.5                  
##  [17] tzdb_0.4.0                  httpuv_1.6.14              
##  [19] assertthat_0.2.1            TCseq_1.22.6               
##  [21] xfun_0.42                   hms_1.1.3                  
##  [23] jquerylib_0.1.4             evaluate_0.23              
##  [25] promises_1.2.1              fansi_1.0.6                
##  [27] caTools_1.18.2              htmlwidgets_1.6.4          
##  [29] DBI_1.2.2                   geneplotter_1.76.0         
##  [31] ellipsis_0.3.2              RSpectra_0.16-1            
##  [33] QFeatures_1.8.0             backports_1.4.1            
##  [35] fontLiberation_0.1.0        prismatic_1.1.1            
##  [37] annotate_1.76.0             fontBitstreamVera_0.1.1    
##  [39] vctrs_0.6.5                 quantreg_5.97              
##  [41] abind_1.4-5                 cachem_1.0.8               
##  [43] withr_3.0.0                 itertools_0.1-3            
##  [45] GenomicAlignments_1.34.1    fdrtool_1.2.17             
##  [47] MultiAssayExperiment_1.24.0 mnormt_2.1.1               
##  [49] cluster_2.1.6               lazyeval_0.2.2             
##  [51] crayon_1.5.2                crul_1.4.0                 
##  [53] labeling_0.4.3              glmnet_4.1-8               
##  [55] edgeR_3.40.2                pkgconfig_2.0.3            
##  [57] nlme_3.1-164                rlang_1.1.3                
##  [59] lifecycle_1.0.4             miniUI_0.1.1.1             
##  [61] MatrixModels_0.5-3          downloader_0.4             
##  [63] fontquiver_0.2.1            httpcode_0.3.0             
##  [65] affyio_1.68.0               extrafontdb_1.0            
##  [67] randomForest_4.7-1.1        rngtools_1.5.2             
##  [69] Matrix_1.6-5                carData_3.0-5              
##  [71] GlobalOptions_0.1.2         png_0.1-8                  
##  [73] rjson_0.2.21                bitops_1.0-7               
##  [75] KernSmooth_2.23-22          Biostrings_2.66.0          
##  [77] blob_1.2.4                  doRNG_1.8.6                
##  [79] shape_1.4.6.1               ggsignif_0.6.4             
##  [81] scales_1.3.0                memoise_2.0.1              
##  [83] magrittr_2.0.3              plyr_1.8.9                 
##  [85] zlibbioc_1.44.0             compiler_4.2.1             
##  [87] pcaMethods_1.90.0           clue_0.3-65                
##  [89] Rsamtools_2.14.0            cli_3.6.2                  
##  [91] affy_1.76.0                 XVector_0.38.0             
##  [93] MASS_7.3-60.0.1             tidyselect_1.2.1           
##  [95] vsn_3.66.0                  stringi_1.8.3              
##  [97] highr_0.10                  yaml_2.3.8                 
##  [99] askpass_1.2.0               locfit_1.5-9.9             
## [101] MALDIquant_1.22.2           ggrepel_0.9.5              
## [103] grid_4.2.1                  sass_0.4.9                 
## [105] ggstats_0.5.1               polynom_1.4-1              
## [107] tools_4.2.1                 timechange_0.3.0           
## [109] parallel_4.2.1              circlize_0.4.16            
## [111] rstudioapi_0.15.0           foreach_1.5.2              
## [113] gridExtra_2.3               farver_2.1.1               
## [115] mzID_1.36.0                 Rtsne_0.17                 
## [117] digest_0.6.35               BiocManager_1.30.22        
## [119] shiny_1.8.0                 gfonts_0.2.0               
## [121] car_3.1-2                   broom_1.0.5                
## [123] later_1.3.2                 ncdf4_1.22                 
## [125] httr_1.4.7                  gdtools_0.3.5              
## [127] AnnotationDbi_1.60.2        ComplexHeatmap_2.14.0      
## [129] colorspace_2.1-0            XML_3.99-0.16.1            
## [131] reticulate_1.35.0           umap_0.2.10.0              
## [133] splines_4.2.1               rematch2_2.1.2             
## [135] plotly_4.10.4               systemfonts_1.0.5          
## [137] xtable_1.8-4                jsonlite_1.8.8             
## [139] pillar_1.9.0                htmltools_0.5.7            
## [141] mime_0.12                   glue_1.7.0                 
## [143] fastmap_1.1.1               BiocParallel_1.32.6        
## [145] class_7.3-22                codetools_0.2-19           
## [147] utf8_1.2.4                  lattice_0.22-5             
## [149] bslib_0.6.1                 curl_5.2.1                 
## [151] gtools_3.9.5                openssl_2.1.1              
## [153] Rttf2pt1_1.3.12             survival_3.5-8             
## [155] rmarkdown_2.26              munsell_0.5.0              
## [157] e1071_1.7-14                GetoptLong_1.0.5           
## [159] GenomeInfoDbData_1.2.9      iterators_1.0.14           
## [161] impute_1.72.3               reshape2_1.4.4             
## [163] gtable_0.3.4                extrafont_0.19